/*
 * MultirootSolver.h
 *
 *  Created on: 23 Jun 2011
 *      Author: allan
 */

#ifndef MULTIROOTSOLVER_H_
#define MULTIROOTSOLVER_H_

// Eigen includes
#include <Eigen/Dense>
using namespace Eigen;

// GSL includes
#include <gsl/gsl_multiroots.h>

class MultirootProblem
{
public:
	virtual void Function(const VectorXd& x, VectorXd& F) = 0;
	
	virtual void Jacobian(const VectorXd& x, MatrixXd& J) = 0;
	
	virtual void FunctionJacobian(const VectorXd& x, VectorXd& F, MatrixXd& J) { Function(x, F); Jacobian(x, J); }
};

class MultirootSolver
{
public:
	enum Algorithm
	{                   /// The list of supported algorithms (for more information, check: http://www.gnu.org/software/gsl/manual/html_node/Algorithms-using-Derivatives.html)
		Hybridj,        /// --- GSL name: (gsl_multiroot_fdfsolver_hybridj)  --- ///
		Hybridsj,       /// --- GSL name: (gsl_multiroot_fdfsolver_hybridsj) --- ///
		StandardNewton, /// --- GSL name: (gsl_multiroot_fdfsolver_newton)   --- ///
		ModifiedNewton  /// --- GSL name: (gsl_multiroot_fdfsolver_gnewton)  --- ///
	};
	
public:
	MultirootSolver();
	
	MultirootSolver(unsigned dimension);
	
	~MultirootSolver();
	
	void
	Initialiaze(unsigned dimension);
	
	void
	SetAlgorithm(Algorithm algorithm);
	
	const unsigned
	GetInteration() const;
	
	void
	SetSolverOptions(double tolerance, unsigned maxiterations);
	
	bool
	Iterate(MultirootProblem& problem, VectorXd& x);
	
	bool
	Solve(MultirootProblem& problem, VectorXd& x);
	
public:
	struct Params
	{
		MultirootProblem* problem;
		VectorXd x;
		VectorXd F;
		MatrixXd J;
	};
	
private:
	unsigned dimension;
	
	Params params;
	
	double tolerance;

	unsigned maxiterations;
	
	/// The number of iterations that the solver used to solve the system of non-linear equations (to be restarted every call to the method Solve)
	unsigned iteration;
	
	const gsl_multiroot_fdfsolver_type* type;
	
	gsl_multiroot_fdfsolver* solver;
	
	gsl_multiroot_function_fdf func;
	
	gsl_vector* xGSL;
};

#endif /* MULTIROOTSOLVER_H_ */
